---
title: "Military Experiment Results"
date: "2025-09-22"
output:
  html_document: default
  pdf_document: default
---

In what follows we fully detail the statistical analysis we conducted. 


## Standards of inference
Unless otherwise stated, we employ the following standards of inference: 

* *strong support* if a coefficient of interest (or difference between coefficients) is in the correct direction and significant with 99% confidence, 
* *support* coefficient of interest (or difference between coefficients) is in the correct direction and significant with 95% confidence,
* *moderate support* coefficient of interest (or difference between coefficients) is in the correct direction and significant with 90% confidence 


We intend to plot results with 95% confidence intervals, unless those results would be deeply misleading.


## Packages used:

```{r, echo = T, results = 'hide', error=FALSE, warning=FALSE, message=FALSE}
rm(list = ls())
library(tidyverse)
library(ggplot2)
library(ggrepel)
library(ggpubr)
library(estimatr)
library(dotwhisker)
library(broom)
library(dplyr)
library(sjPlot)
library(ivreg)
library(readxl)
library(gtools)
library(stargazer)
library(cobalt)
library(summarytools)
library(pwrss)
library(pwr)
```

## Importing, Cleaning, & Stacking Data. 

```{r,echo = T, results = 'hide',  error=FALSE, warning=FALSE, message=FALSE}
milexp1 <- read_xlsx("milexp1.xlsx")
milexp2 <- read_xlsx("milexp2.xlsx")
```

```{r,echo = T, results = 'hide',  error=FALSE, warning=FALSE, message=FALSE}
milexp1$T1 = rep(NA, nrow(milexp1))
milexp1$T1[milexp1$group == 1] <- 1 # partisan agreement
milexp1$T1[milexp1$group == 2] <- 1 # partisan agreement
milexp1$T1[milexp1$group == 3] <- 0 # partisan disagreement
milexp1$T1[milexp1$group == 4] <- 0 # partisan disagreement

milexp1$T2 = rep(NA, nrow(milexp1))
milexp1$T2[milexp1$group == 1] <- 1 # incumbent party identifier
milexp1$T2[milexp1$group == 2] <- 0 # opposition party identifier
milexp1$T2[milexp1$group == 3] <- 1 # incumbent party identifier
milexp1$T2[milexp1$group == 4] <- 0 # opposition party identifier

milexp1$female[milexp1$female == 1] <- 2
milexp1$female[milexp1$female == 0] <- 1

milexp1$priorexp[milexp1$priorexp == 0] <- 2

milexp1$study = rep(NA, nrow(milexp1))
milexp1$study <- 1
milexp1$ID <- (1:nrow(milexp1))

milexp1 <- rename(milexp1, "risklife" = "SLE1-risklife")
milexp1 <- rename(milexp1, "risksubordin" = "SLE2-risksubordin")
milexp1 <- rename(milexp1, "discretion" = "SLE3-discretion")
milexp1 <- rename(milexp1, "interunit" = "SLE4-interunit")
milexp1 <- rename(milexp1, "unitcohesion" = "SLE5-unitcohesion")
milexp1 <- rename(milexp1, "obey" = "SLE6-obey")

milexp1 <- rename(milexp1, "AM_1" = "AM1-threat")
milexp1 <- rename(milexp1, "AM_2" = "AM2-partytype")
milexp1 <- rename(milexp1, "AM_3" = "AM3-artificial")


milexp2$T1 = rep(NA, nrow(milexp2))
milexp2$T1[milexp2$group == 1] <- 1 # partisan agreement
milexp2$T1[milexp2$group == 2] <- 1 # partisan agreement
milexp2$T1[milexp2$group == 3] <- 0 # partisan disagreement
milexp2$T1[milexp2$group == 4] <- 0 # partisan disagreement

milexp2$T2 = rep(NA, nrow(milexp2))
milexp2$T2[milexp2$group == 1] <- 1 # incumbent party identifier
milexp2$T2[milexp2$group == 2] <- 0 # opposition party identifier
milexp2$T2[milexp2$group == 3] <- 1 # incumbent party identifier
milexp2$T2[milexp2$group == 4] <- 0 # opposition party identifier

milexp2$female[milexp2$female == 1] <- 2
milexp2$female[milexp2$female == 0] <- 1

milexp2$study = rep(NA, nrow(milexp2))
milexp2$study <- 2
milexp2$ID <- ((nrow(milexp1)+1):(nrow(milexp2)+nrow(milexp1)))

milexp2 <- rename(milexp2, "risklife" = "SLE1-risklife")
milexp2 <- rename(milexp2, "risksubordin" = "SLE2-risksubordin")
milexp2 <- rename(milexp2, "discretion" = "SLE3-discretion")
milexp2 <- rename(milexp2, "interunit" = "SLE4-interunit")
milexp2 <- rename(milexp2, "unitcohesion" = "SLE5-unitcohesion")
milexp2 <- rename(milexp2, "obey" = "SLE6-obey")

milexp2 <- rename(milexp2, "AM_1" = "AM1-threat")
milexp2 <- rename(milexp2, "AM_2" = "AM2-partytype")
milexp2 <- rename(milexp2, "AM_3" = "AM3-artificial")
```

```{r, error=FALSE, warning=FALSE, message=FALSE}
milexp1.1 <- milexp1 %>% 
  select(risklife:obey, NIM:partyIDstrength, anger:ID, group) %>%
  group_by(ID)
milexp2.1 <- milexp2 %>% 
  select(NIM:AM_3, age:partyIDstrength, anger:group, T1:ID) %>%
  group_by(ID)

milexp <- rbind(milexp1.1, milexp2.1)

datstack <- milexp %>%
  pivot_longer(cols = "risklife":"obey", 
               names_to = "Variable", 
               values_to = "Value",
               values_drop_na = F)
```


```{r, error=FALSE, warning=FALSE, message=FALSE}
print.data.frame(head(datstack))
``` 

Note here the `value` column represents the outcomes, and the `Variable` column is the label of the particular outcome.   

The key feature of this data structure is that our predictions point in the same direction for each outcome. That is, we typically expect `T1` to predict higher number for all outcomes, and a stronger effect for all outcomes when `T2=1`.
  
  
# Hypothesis 1 and 2

In sections 3.4 and 3.6 of the manuscript, we describe 6 tests across our main two hypotheses. For each hypothesis, we examine one *global test* that examines the results across both study contexts and all 6 treatment outcomes. Then two *context-specific tests* that isolate a specific context (study 1 is the North Korean collapse Scenario, study 2 is the US alliance scenario). 

## Hypothesis 1

The coefficient of interest is T1 (partisan (dis)agreement) and we use our main standard of inference.

Our global test
```{r, error=FALSE, warning=FALSE, message=FALSE}
H1 <- lm_robust(Value ~ T1 + Variable + study,  
                clusters = ID,
                data = datstack)

summary(H1)

``` 


Our context-specific test
```{r, error=FALSE, warning=FALSE, message=FALSE}

H11 <-   datstack %>% filter(study==1) %>% 
  lm_robust(Value ~ T1 + Variable,  
                clusters = ID,
                data = .)
summary(H11)


H12 <- datstack %>% filter(study==2) %>% 
  lm_robust(Value ~ T1 + Variable,  
                clusters = ID,
                data = .)

summary(H12)

```

### A note on conflicting results
If we find support for the global test at 95% and at least moderate support for both context-specific test, we say that we have support for H1 across contexts. We find moderate support (significant at the .90 level) for the main hypothesis.

It is possible that the global and context-specific tests conflict. If we find no support for one local test and support for the other we say that we find support for H1 in one particular context but not the other. 

**Note:** We apply the same reasoning for H1 and H2. 

## H2: Interaction of treatments

```{r, error=FALSE, warning=FALSE, message=FALSE}
## Global Model

H2 <- lm_robust(Value ~ T1*T2 + Variable + study,
                clusters = ID,
                data = datstack)
summary(H2)

## Experiment 1

H21 <-  datstack %>% filter(study==1) %>%  
  lm_robust(Value ~ T1*T2 + Variable,
                clusters = ID,
                data = .)

## Experiment 2

H22 <-   datstack %>% filter(study==2) %>%  
  lm_robust(Value ~ T1*T2 + Variable,
                clusters = ID,
                data = .)

```

We first tried to interpret the results from a visual inspection of the following marginal effects plots, as stated in our pre-analysis plan.

```{r, error=FALSE, warning=FALSE, message=FALSE}

## Global Model
datstack$Elite <- rep(NA, nrow(datstack))
datstack$Elite[datstack$T1 == 0] <- "Disagree"
datstack$Elite[datstack$T1 == 1] <- "Agree"
datstack$Partisanship[datstack$T2 == 0] <- "Opposition"
datstack$Partisanship[datstack$T2 == 1] <- "Incumbent"

H2.0 <- lm_robust(Value ~ Elite*Partisanship + Variable + study,  
            clusters = ID,
            data = datstack)

summary(H2.0)

## Experiment 1

datA <- datstack %>% filter(study==1)

H2.1 <- lm_robust(Value ~ Elite*Partisanship + Variable,  
            clusters = ID,
            data = datA)

summary(H2.1)
 
## Experiment 2

datB <- datstack %>% filter(study==2)

H2.2 <- lm_robust(Value ~ Elite*Partisanship + Variable,
            clusters = ID,
            data = datB)
summary(H2.2)


###### PLOTS
## Interaction Plots:

# Global Model
# For the plot_model() function, we always use "discretion" as a baseline category
set_theme(base = theme_classic())
p1 <- plot_model(H2.0, type = "pred", terms = c("Partisanship", "Elite"), colors = "bw",
           condition = c(Variable = "discretion"), title = "Global Model", 
           dot.size = 1, axis.title = c("Partisanship", "Military Cohesion")) + 
      theme(text = element_text(size = 6)) + annotate("text", x=2, y=4.2, 
                                                      label = "* p = .042",
                                                      size = 2, color = "red2") +
      annotate("text", x=1.05, y=4.19, label = "* p = .346", size = 2) +
      annotate("text", x=2, y=4.212, label = "diff = .111", size = 2) +
      annotate("text", x=1.05, y=4.202, label = "diff = .051", size = 2)

# Experiment 1
p2 <- plot_model(H2.1, type = "pred", terms = c("Partisanship", "Elite"), colors = "bw", 
           condition = c(Variable = "discretion"), title = "Experiment 1",
           dot.size = 1, axis.title = c("Partisanship", "Military Cohesion")) + 
      theme(text = element_text(size = 6)) + annotate("text", x=2, y=3.515, 
                                                      label = "* p = .339",
                                                      size = 2) +
      annotate("text", x=1.05, y=3.47, label = "* p = .712", size = 2) +
      annotate("text", x=2, y=3.53, label = "diff = .076", size = 2) +
      annotate("text", x=1.05, y=3.485, label = "diff = .029", size = 2)

# Experiment 2
p3 <- plot_model(H2.2, type = "pred", terms = c("Partisanship", "Elite"), colors = "bw", 
           condition = c(Variable = "discretion"), title = "Experiment 2",
           dot.size = 1, axis.title = c("Partisanship", "Military Cohesion")) + 
      theme(text = element_text(size = 6)) + annotate("text", x=2, y=4.905, 
                                                      label = "* p = .072",
                                                      size = 2, color = "red2") +
      annotate("text", x=1.05, y=4.925, label = "* p = .694", size = 2) +
      annotate("text", x=2, y=4.92, label = "diff = .118", size = 2) +
      annotate("text", x=1.05, y=4.94, label = "diff = .025", size = 2)

ggarrange(p1, p2, p3, ncol = 3, nrow = 1, common.legend = T)

```

As there is a small amount of overlap such that we cannot draw an inference from visual inspection, we decided to proceed with the standard analytical tests (see below; not preregistered tests).

# Dimension-specific analysis. 

We now present the plots from the separate analysis of experiments 1 and 2, along with the global model result.

```{r, error=FALSE, warning=FALSE, message=FALSE}

names <- c("risklife", "risksubordin", "discretion", "interunit", "unitcohesion",
           "obey")

datstack1 <- datstack[!is.na(datstack$T1), ]


H11n_t <- data.frame()
for(i in 1:6){
  H11n_t <- rbind(H11n_t, datstack1 %>% filter(Variable == names[i]) %>% 
  lm(Value ~ T1*T2, data = .) %>% 
  tidy() %>% filter(term== "T1") %>% 
  mutate(model = names[i]))
}


H11n_t1 <- data.frame()
for(i in 1:6){
  H11n_t1 <- rbind(H11n_t1, datstack1 %>% filter(study==1 & Variable == names[i]) %>% 
  lm(Value ~ T1*T2, data = .) %>% 
  tidy() %>% filter(term== "T1") %>% 
  mutate(model = names[i]))
}


H11n_t2 <- data.frame()
for(i in 1:6){
  H11n_t2 <- rbind(H11n_t2, datstack1 %>% filter(study==2 & Variable == names[i]) %>% 
  lm(Value ~ T1*T2, data = .) %>% 
  tidy() %>% filter(term== "T1") %>% 
  mutate(model = names[i]))
}


H11n_t
H11n_t1
H11n_t2



p_t <- ggplot(H11n_t,
               aes(x = estimate,
                   y = model,
                   color = model,
                   shape = model)) +
  geom_vline(xintercept = 0, linetype = 2, color = "black") +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - 1.96 * std.error,
                     xmax = estimate + 1.96 * std.error),
                 height = 0.1) +
  theme_bw() +
  xlab("Global Model") +
  ylab("") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  scale_color_grey(
    start = .1, end = .1,
    name = "Military Cohesion",
    labels = c("Obeyorders", "Unitcohesion", "Interunit",
               "Discretion", "Risksubordin", "Risklife")
  ) +
  scale_shape_manual(
    name = "Military Cohesion",
    values = 1:6,
    labels = c("Obeyorders", "Unitcohesion", "Interunit",
               "Discretion", "Risksubordin", "Risklife")
  ) +
  guides(color = guide_legend(order = 1),
         shape = guide_legend(order = 1))


p_t1 <- ggplot(H11n_t1,
               aes(x = estimate,
                   y = model,
                   color = model,
                   shape = model)) +
  geom_vline(xintercept = 0, linetype = 2, color = "black") +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - 1.96 * std.error,
                     xmax = estimate + 1.96 * std.error),
                 height = 0.1) +
  theme_bw() +
  xlab("Experiment 1") +
  ylab("") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  scale_color_grey(
    start = .1, end = .1,
    name = "Military Cohesion",
    labels = c("Obeyorders", "Unitcohesion", "Interunit",
               "Discretion", "Risksubordin", "Risklife")
  ) +
  scale_shape_manual(
    name = "Military Cohesion",
    values = 1:6,
    labels = c("Obeyorders", "Unitcohesion", "Interunit",
               "Discretion", "Risksubordin", "Risklife")
  ) +
  guides(color = guide_legend(order = 1),
         shape = guide_legend(order = 1))


p_t2 <- ggplot(H11n_t2,
               aes(x = estimate,
                   y = model,
                   color = model,
                   shape = model)) +
  geom_vline(xintercept = 0, linetype = 2, color = "black") +
  geom_point(size = 3) +
  geom_errorbarh(aes(xmin = estimate - 1.96 * std.error,
                     xmax = estimate + 1.96 * std.error),
                 height = 0.1) +
  theme_bw() +
  xlab("Experiment 2") +
  ylab("") +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  scale_color_grey(
    start = .1, end = .1,
    name = "Military Cohesion",
    labels = c("Obeyorders", "Unitcohesion", "Interunit",
               "Discretion", "Risksubordin", "Risklife")
  ) +
  scale_shape_manual(
    name = "Military Cohesion",
    values = 1:6,
    labels = c("Obeyorders", "Unitcohesion", "Interunit",
               "Discretion", "Risksubordin", "Risklife")
  ) +
  guides(color = guide_legend(order = 1),
         shape = guide_legend(order = 1))


ggarrange(p_t, p_t1, p_t2, ncol = 3, nrow = 1, common.legend = T)

```

# Heterogenous effect analysis

We measure partisan variables, military experience variables in both studies. But we only measure the full battery of biographical attributes in study 1, and we focus our analysis here. 

## Data manipulations
Overall, we collect data on 11 psychological and demographic attributes. In the manuscript, we intend to plot the moderated treatment effects of a few, while controlling for the rest. We now explain how we manipulate the data for the variables of interest.

The strength of party identification uses the same 5-point scale used in the survey questionnaire.

Region variable, which represents the province in which a respondent's hometown is located, is a dichotomous variable, which has value 2 and 3 for two provinces with historically strong partisan animosity against each other, Gyeongsang and Jeolla provinces, respectively, and 1 for the others.

Both Gender and School Year are dichotomous. Gender variable has value 1 if a respondent is male and 2 if a respondent is female. School Year variable similarly has value 1 if a respondent is a junior and 2 if a respondent is a senior.

Psychological variables of particular interest, in relation to previous works, are Anger and Risk Propensity. Both are measured by 7-point Likert scale in the survey and will be transformed to be ordinal variables with three levels. The levels of Anger are divided into low (values 1, 2, 3), moderate (value 4), and high (values 5, 6, 7) and the levels of Risk Propensity into risk-verse (values 1, 2, 3), risk-neutral (value 4), and risk-taking (values 5, 6, 7). The smaller the value in the agree/disagree 7-point Likert scale question is, the more likely a respondent agrees to the given statement in the question.

Other variables included for control are the level of education (1 if high school, 2 if Associate's, 3 if Bachelor's), prior military experience (1 if no, 2 if yes), the number of siblings (1 if none, 2 if one, 3 if two, 4 if three or more), moral sense or consciousness (low (values 1, 2, 3), neutral (value 4), high (values 5, 6, 7)), military assertiveness (dove (values 1, 2, 3), neutral (value 4), hawk (values 5, 6, 7)) and age (1 if below 25 and 2 if above or equal to 25).

* In each of the heterogeneous treatment effect models outlined below, the moderator of the model's interest will be operationalized with this transformed factor variable version to make interpretation easy. Other moderators that are used as control variables in the model will be treated as continuous variables unless they are categorical by nature (e.g. region, gender). This is mainly to not decrease too much statistical power by unnecessarily adding too many factor variables.

```{r, error=FALSE, warning=FALSE, message=FALSE}
###############
#   Heterogeneous Treatment Effects: Moderator Manipulation
##################

## Experiment 1

datA$Elite <- rep(NA, nrow(datA))
datA$Party_ID_Strength <- rep(NA, nrow(datA))
datA$Region <- rep(NA, nrow(datA))
datA$Female <- rep(NA, nrow(datA))
datA$School_Yr <- rep(NA, nrow(datA))
datA$Anger <- rep(NA, nrow(datA))
datA$Risk_Propensity <- rep(NA, nrow(datA))

datA$Elite[datA$T1 == 1] <- "Partisan Agreement"
datA$Elite[datA$T1 == 0] <- "Partisan Disagreement"

datA$Party_ID_Strength[datA$partyIDstrength == 1] <- "Independent"
datA$Party_ID_Strength[datA$partyIDstrength == 2] <- "Weak"
datA$Party_ID_Strength[datA$partyIDstrength == 3] <- "Moderate"
datA$Party_ID_Strength[datA$partyIDstrength == 4] <- "Strong"
datA$Party_ID_Strength[datA$partyIDstrength == 5] <- "Very Strong"
datA$Party_ID_Strength <- factor(datA$Party_ID_Strength,
                                 levels=c("Independent", "Weak", "Moderate",
                                          "Strong", "Very Strong"))

datA$Region[datA$hometown == 1] <- "Others"
datA$Region[datA$hometown == 2] <- "Gyeongsang"
datA$Region[datA$hometown == 3] <- "Jeolla"

datA$Female[datA$female == 1] <- "Male"
datA$Female[datA$female == 2] <- "Female"

datA$School_Yr[datA$rank == 1] <- "Junior"
datA$School_Yr[datA$rank == 2] <- "Senior"

datA$Anger[datA$anger < 4] <- "Low"
datA$Anger[datA$anger == 4] <- "Moderate"
datA$Anger[datA$anger > 4] <- "High"
datA$Anger <- factor(datA$Anger, levels=c("Low", "Moderate", "High"))

datA$Risk_Propensity[datA$risk < 4] <- "Risk-Averse"
datA$Risk_Propensity[datA$risk == 4] <- "Risk-Neutral"
datA$Risk_Propensity[datA$risk > 4] <- "Risk-Taking"
datA$Risk_Propensity <- factor(datA$Risk_Propensity, levels=c("Risk-Averse",
                                                              "Risk-Neutral", 
                                                              "Risk-Taking"))

datA$Education <- rep(NA, nrow(datA))
datA$Siblings <- rep(NA, nrow(datA))
datA$Moral <- rep(NA, nrow(datA))
datA$Dove <- rep(NA, nrow(datA))
datA$Age <- rep(NA, nrow(datA))

datA$Education[datA$edu == 1] <- "High School"
datA$Education[datA$edu == 2] <- "A.A."
datA$Education[datA$edu == 3] <- "B.A."

datA$Siblings[datA$sibling == 1] <- "0"
datA$Siblings[datA$sibling == 2] <- "1"
datA$Siblings[datA$sibling == 3] <- "2"
datA$Siblings[datA$sibling == 4] <- "3 or more"

datA$Moral[datA$moral < 4] <- "Low"
datA$Moral[datA$moral == 4] <- "Neutral"
datA$Moral[datA$moral > 4] <- "High"
datA$Moral <- factor(datA$Moral, levels=c("Low", "Neutral", "High"))

datA$Dove[datA$dove < 4] <- "Hawk"
datA$Dove[datA$dove == 4] <- "Neutral"
datA$Dove[datA$dove > 4] <- "Dove"
datA$Dove <- factor(datA$Dove, levels=c("Dove", "Neutral", "Hawk"))

datA$Age[datA$age < 25] <- "Below 25"
datA$Age[datA$age >= 25] <- "Above or Equal to 25"


## Experiment 2

datC <- milexp2 %>% 
  select(NIM:AM_3, age:partyIDstrength, anger:group, T1:ID, class) %>%
  group_by(ID)

datCstack <- datC %>%
  pivot_longer(cols = "risklife":"obey", 
               names_to = "Variable", 
               values_to = "Value",
               values_drop_na = F)

datCstack$Elite <- rep(NA, nrow(datCstack))
datCstack$Elite[datCstack$T1 == 1] <- "Partisan Agreement"
datCstack$Elite[datCstack$T1 == 0] <- "Partisan Disagreement"

datCstack$Class <- rep(NA, nrow(datCstack))
datCstack$Class[datCstack$class == "officer1"] <- "Active Duty Officer"
datCstack$Class[datCstack$class == "officer2"] <- "Active Duty Officer"
datCstack$Class[datCstack$class == "conscript"] <- "Conscripted Soldier"
datCstack$Class[datCstack$class == "cadet"] <- "Cadet"
```

We run stacked regressions for estimating heterogeneous treatment effects and plot marginal effects as follows.

```{r, error=FALSE, warning=FALSE, message=FALSE}
## Party ID Strength
hte11 <- lm_robust(Value ~ Elite*Party_ID_Strength + Region + Female + School_Yr 
                   + edu + sibling + age + Anger + Risk_Propensity 
                   + moral + dove + Variable,  
            clusters = ID,
            data = datA)

set_theme(base = theme_classic())
p11<- plot_model(hte11, type = "pred", terms = c("Party_ID_Strength", "Elite"), 
           colors = "bw", condition = c(Variable = "discretion"),
           title = "Party ID Strength (Exp1)", dot.size = 1,
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5))



## Region
hte12 <- lm_robust(Value ~ Elite*Region + Party_ID_Strength + Female + School_Yr 
                  + edu + sibling + age + Anger + Risk_Propensity 
                  + moral + dove + Variable,  
            clusters = ID,
            data = datA)

p12 <- plot_model(hte12, type = "pred", terms = c("Region", "Elite"), colors = "bw", 
           title = "Region (Exp1)", dot.size = 1, condition = c(Variable = "discretion"),
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5))



## Gender
hte13 <- lm_robust(Value ~ Elite*Female + Party_ID_Strength + Region + School_Yr 
                  + edu + sibling + age + Anger + Risk_Propensity 
                   + moral + dove + Variable,  
            clusters = ID,
            data = datA)

p13 <- plot_model(hte13, type = "pred", terms = c("Female", "Elite"), colors = "bw", 
           title = "Gender (Exp1)", dot.size = 1, condition = c(Variable = "discretion"),
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5))


## School Year
hte14 <- lm_robust(Value ~ Elite*School_Yr + Party_ID_Strength + Region + Female 
                   + edu + sibling + age + Anger + Risk_Propensity 
                   + moral + dove + Variable,  
            clusters = ID,
            data = datA)

p14 <- plot_model(hte14, type = "pred", terms = c("School_Yr", "Elite"), colors = "bw", 
           title = "School Year (Exp1)", dot.size = 1,
           condition = c(Variable = "discretion"),
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5)) +
       annotate("text", x=1.05, y=3.5, label = "* p = .118", size = 2)


## Anger
hte15 <- lm_robust(Value ~ Elite*Anger + Party_ID_Strength + Region + Female 
                   + School_Yr + edu + sibling + age + Risk_Propensity 
                   + moral + dove + Variable,  
            clusters = ID,
            data = datA)

p15 <- plot_model(hte15, type = "pred", terms = c("Anger", "Elite"), colors = "bw", 
           title = "Anger (Exp1)", dot.size = 1, condition = c(Variable = "discretion"),
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5))



## Risk Propensity
hte16 <- lm_robust(Value ~ Elite*Risk_Propensity + Party_ID_Strength + Region 
                   + Female + School_Yr + edu + sibling + age + Anger 
                   + moral + dove + Variable,  
            clusters = ID,
            data = datA)

p16 <- plot_model(hte16, type = "pred", terms = c("Risk_Propensity", "Elite"), 
           colors = "bw", title = "Risk Propensity (Exp1)", dot.size = 1,
           condition = c(Variable = "discretion"),
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5))


## Morality
hte17 <- lm_robust(Value ~ Elite*Moral + Party_ID_Strength + Region 
                   + Female + School_Yr + edu + sibling + age + Anger 
                   + Risk_Propensity + dove + Variable,  
            clusters = ID,
            data = datA)

p17 <- plot_model(hte17, type = "pred", terms = c("Moral", "Elite"), colors = "bw", 
           title = "Moral Sense (Exp1)", dot.size = 1,
           condition = c(Variable = "discretion"),
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5))


## Militant Assertiveness
hte18 <- lm_robust(Value ~ Elite*Dove + Party_ID_Strength + Region 
                   + Female + School_Yr + edu + sibling + age + Anger 
                   + moral + Risk_Propensity + Variable,  
            clusters = ID,
            data = datA)

p18 <- plot_model(hte18, type = "pred", terms = c("Dove", "Elite"), colors = "bw", 
           title = "Militant Assertiveness (Exp1)", dot.size = 1,
           condition = c(Variable = "discretion"),
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5)) +
       annotate("text", x=1.05, y=3.5, label = "* p = .112", size = 2)


## Military Class
hte21 <- lm_robust(Value ~ Elite*Class + age + female + rank + edu + hometown +
                   sibling + partyIDstrength + anger + moral + risk + dove +
                   Variable,  
            clusters = ID,
            data = datCstack)

p21 <- plot_model(hte21, type = "pred", terms = c("Class", "Elite"), colors = "bw", 
           title = "Military Class (Exp2)", dot.size = 1,
           condition = c(Variable = "discretion"),
           axis.title = c("", "Military Cohesion")) + theme(text = element_text(size = 5))


ggarrange(p11, p12, p13, p14, p15, p16, p17, p18, p21, ncol = 3, nrow = 3,
          common.legend = T)

```

Testing the mechanisms

```{r, error=FALSE, warning=FALSE, message=FALSE}
####################### 
#  Testing the mechanisms. 
####################

#### Direct tests:
# As a first cut, we make sure that our treatments impact our mechanisms 
# by estimating their effects on the questions that most closely reflect the mechanisms. 

## Here is the direct test (H1): 
# From treatment to mediators
NIM <- lm(NIM ~ T1, milexp)
TIM <- lm(TIM ~ T1, milexp)

summary(NIM)
summary(TIM)


# From mediators to DV
NIM2 <- lm_robust(Value ~ NIM + Variable + study,  
                clusters = ID,
                data = datstack)
TIM2 <- lm_robust(Value ~ TIM + Variable + study,  
                clusters = ID,
                data = datstack)

summary(NIM2)
summary(TIM2)


## If these are positive and significant, we may proceed to more complex testing 
# and visual inspection if we also pass other benchmarks. 

## If they are not, we note that we did not find support for the impact of treatment 
# on a particular mechanism. 


## If we ALSO observe a difference for the effect of T1|T2, then we may proceed to 
# more complex testing and visual inspection if we also pass other benchmarks. 

#If they are not, we note that we did not find support for the impact of treatment 
# on a particular mechanism. 

IV1 <- ivreg(Value ~ NIM + TIM + Variable + study | T1 + T2 + T1:T2 + Variable,
                clusters = ID,
                data = datstack)

summary(IV1)

## Direct test: 
M1 <- lm_robust(Value ~ NIM + TIM + Variable + study,
                clusters = ID,
                data = datstack)

summary(M1)


## Alternative Mechanisms
AM1 <- lm(AM_1 ~T1, milexp)
AM2 <- lm(AM_2 ~T1, milexp)
AM3 <- lm(AM_3 ~T1, milexp)

summary(AM1)
summary(AM2)
summary(AM3)


AM <- lm_robust(Value ~ AM_1 + AM_2 + AM_3 + Variable + study,
                clusters = ID,
                data = datstack)

summary(AM)
```

From here, we include complementary statistical tests that we did not preregister. These include:

(1) Analytical tests for H2: We noted in our preregistered Rmarkdown file that we would conduct analytical tests for interaction effects if we cannot clearly identify the effects with the marginal plot. As this has turned out to be the case for the opposition party identifiers, we report results from the one-sided t-test, which is commonly used as an analytical test comparing marginal effects. Opposition group identifiers show statistically significant difference in military cohesion depending on the partisan disagreement treatment (at the .90 level). We annotate p-values for the significant results in our marginal plot.

```{r, error=FALSE, warning=FALSE, message=FALSE}
###### ANALYTICAL TESTS (Two-Tailed T-Tests)
## One-Sided Test (Incumbent Side)
# Global Model
h2full.1 <- t.test(datstack$Value[datstack$T2 == 1 & datstack$T1 == 1],
       datstack$Value[datstack$T2 == 1 & datstack$T1 == 0],
       alternative = "two.sided", conf.level = .90)
h2full.1

# Experiment 1
h2st1.1 <- t.test(datA$Value[datA$T2 == 1 & datA$T1 == 1],
       datA$Value[datA$T2 == 1 & datA$T1 == 0],
       alternative = "two.sided", conf.level = .90)
h2st1.1

# Experiment 2
h2st2.1 <- t.test(datB$Value[datB$T2 == 1 & datB$T1 == 1],
       datB$Value[datB$T2 == 1 & datB$T1 == 0],
       alternative = "two.sided", conf.level = .90)
h2st2.1


## One-Sided Test (Opposition Side)
# Global Model: Significant at the .90 level
h2full.2 <- t.test(datstack$Value[datstack$T2 == 0 & datstack$T1 == 1],
       datstack$Value[datstack$T2 == 0 & datstack$T1 == 0],
       alternative = "two.sided", conf.level = .90)
h2full.2 

# Experiment 1
h2st1.2 <- t.test(datA$Value[datA$T2 == 0 & datA$T1 == 1],
       datA$Value[datA$T2 == 0 & datA$T1 == 0],
       alternative = "two.sided", conf.level = .90)
h2st1.2

# Experiment 2: Significant at the .90 level
h2st2.2 <- t.test(datB$Value[datB$T2 == 0 & datB$T1 == 1],
       datB$Value[datB$T2 == 0 & datB$T1 == 0],
       alternative = "two.sided", conf.level = .90)
h2st2.2
```

(2) Analytical tests for heterogeneous treatment effects (HTE): We report results from the one-sided t-tests for HTE that cannot be clearly judged by exploring the marginal plots. These include the gender and school year variables, and risk-averse,
low moral sense, and dove groups. Junior cadets and doves in partisan agreement group and those in partisan disagreement group are found to show difference in military cohesion that is significant at the .90 level. We annotate p-values for the significant results in our marginal plot.

```{r, error=FALSE, warning=FALSE, message=FALSE}
### Analytical Tests (Two-Tailed T-Tests) for HTE
## Gender
# Male - insignificant
t.test(datA$Value[datA$Female == "Male" & datA$T1 == 1], 
       datA$Value[datA$Female == "Male" & datA$T1 == 0], alternative = "two.sided",
       conf.level = .90)

# Female - insignificant
t.test(datA$Value[datA$Female == "Female" & datA$T1 == 1], 
       datA$Value[datA$Female == "Female" & datA$T1 == 0], alternative = "two.sided",
       conf.level = .90)


## School Year
# Junior
t.test(datA$Value[datA$School_Yr == "Junior" & datA$T1 == 1], 
       datA$Value[datA$School_Yr == "Junior" & datA$T1 == 0], alternative = "two.sided",
       conf.level = .90)

# Senior - insignificant
t.test(datA$Value[datA$School_Yr == "Senior" & datA$T1 == 1], 
       datA$Value[datA$School_Yr == "Senior" & datA$T1 == 0], alternative = "two.sided",
       conf.level = .90)


## Risk Propensity (for risk-averse) - insignificant
t.test(datA$Value[datA$Risk_Propensity == "Risk-Averse" & datA$T1 == 1],
       datA$Value[datA$Risk_Propensity == "Risk-Averse" & datA$T1 == 0], 
       alternative = "two.sided", conf.level = .90)


## Moral Sense (for low morality) - insignificant
t.test(datA$Value[datA$Moral == "Low" & datA$T1 == 1], 
       datA$Value[datA$Moral == "Low" & datA$T1 == 0], alternative = "two.sided",
       conf.level = .90)


## Militant Assertiveness (for doves)
t.test(datA$Value[datA$Dove == "Dove" & datA$T1 == 1], 
       datA$Value[datA$Dove == "Dove" & datA$T1 == 0], alternative = "two.sided",
       conf.level = .90)
```

(3) We provide descriptive summary plots for responses to the alternative mechanism (party type plausibility (AM_2), scenario realism (AM_3)) questions. The majority of respondents across all treatment groups (T1, T2) equally believe (i) that the party they favor is plausibly expected to take the position assigned in the scenario they read and (ii) that the scenario presented is sufficiently realistic.

```{r, error=FALSE, warning=FALSE, message=FALSE}
## Descriptive Representation of Party Type Plausibility by Experimental Conditions
milexp$Elite[milexp$T1 == 0] <- "Disagree"
milexp$Elite[milexp$T1 == 1] <- "Agree"
milexp$Partisanship[milexp$T2 == 0] <- "Opposition"
milexp$Partisanship[milexp$T2 == 1] <- "Incumbent"

milexp %>%
  mutate(name1 = fct_relevel(Elite, 
            "Agree", "Disagree")) %>%
  ggplot(mapping = aes(x = AM_2, fill = name1)) +
  geom_bar(mapping = aes(fill = name1, position = "dodge"), alpha = .6) + 
  xlab("Party Type Plausibility") + ylab("n of responses") +
  scale_fill_manual(name = "Elite", values=c("skyblue2", "pink2")) +
  scale_x_continuous(breaks = seq(1,7,1), labels=c("1" = "Very Implausible",
                            "2" = "Implausible", "3" = "Somewhat Implausible", 
                            "4" = "Neither", "5" = "Somewhat Plausible", 
                            "6" = "Plausible", "7" = "Very Plausible")) +
  theme(axis.text.x = element_text(angle = 60, size = 10, hjust = 1, vjust = 1))

milexp %>%
  mutate(name2 = fct_relevel(Partisanship, 
            "Incumbent", "Opposition")) %>%
  ggplot(mapping = aes(x = AM_2, fill = name2)) +
  geom_bar(mapping = aes(fill = name2, position = "dodge"), alpha = .6) + 
  xlab("Party Type Plausibility") + ylab("n of responses") +
  scale_fill_manual(name = "Partisanship", values=c("purple2", "gold2")) +
  scale_x_continuous(breaks = seq(1,7,1), labels=c("1" = "Very Implausible",
                            "2" = "Implausible", "3" = "Somewhat Implausible",
                            "4" = "Neither", "5" = "Somewhat Plausible",
                            "6" = "Plausible", "7" = "Very Plausible")) +
  theme(axis.text.x = element_text(angle = 60, size = 10, hjust = 1, vjust = 1))


## Descriptive Representation of Scenario Realism by Experimental Conditions
milexp %>%
  mutate(name1 = fct_relevel(Elite, 
            "Agree", "Disagree")) %>%
  ggplot(mapping = aes(x = AM_3, fill = name1)) +
  geom_bar(mapping = aes(fill = name1, position = "dodge"), alpha = .6) + 
  xlab("Scenario Realism") + ylab("n of responses") +
  scale_fill_manual(name = "Elite", values=c("skyblue2", "pink2")) +
  scale_x_continuous(breaks = seq(1,7,1), labels=c("1" = "Very Unrealistic",
                            "2" = "Unrealistic", "3" = "Somewhat Unrealistic",
                            "4" = "Neither", "5" = "Somewhat Realistic",
                            "6" = "Realistic", "7" = "Very Realistic")) +
  theme(axis.text.x = element_text(angle = 60, size = 10, hjust = 1, vjust = 1))

milexp %>%
  mutate(name2 = fct_relevel(Partisanship, 
            "Incumbent", "Opposition")) %>%
  ggplot(mapping = aes(x = AM_3, fill = name2)) +
  geom_bar(mapping = aes(fill = name2, position = "dodge"), alpha = .6) + 
  xlab("Scenario Realism") + ylab("n of responses") +
  scale_fill_manual(name = "Partisanship", values=c("purple2", "gold2")) +
  scale_x_continuous(breaks = seq(1,7,1), labels=c("1" = "Very Unrealistic", 
                              "2" = "Unrealistic", "3" = "Somewhat Unrealistic",
                              "4" = "Neither", "5" = "Somewhat Realistic",
                              "6" = "Realistic", "7" = "Very Realistic")) +
  theme(axis.text.x = element_text(angle = 60, size = 10, hjust = 1, vjust = 1))
```

(4) Auxiliary test results: We report descriptive statistics for our global model and each wave of the experiments. We also share balance test results for the global model, where we check the balance of covariates across the main treatment groups (partisan agreement vs. disagreement). The hometown province variable is slightly imbalanced, as its mean standardized difference between the two treatment groups is more than .10 distance away from zero. Yet, in the main analysis of H1 provided earlier, we control for this variable and the result does not change depending on the adjustment.

```{r, error=FALSE, warning=FALSE, message=FALSE}
####################### 
#  Auxiliary Tests
####################

### Descriptive Statistics
## Total
# Global Model
milexp_a <- bind_rows(milexp1, milexp2)
print(dfSummary(milexp_a), method = "render")

# Experiment 1
print(dfSummary(milexp1), method = "render")

# Experiment 2
print(dfSummary(milexp2), method = "render")


## T1 = 1
# Experiment 1
print(dfSummary(milexp1[milexp1$T1 == 1,]), method = "render")

# Experiment 2
print(dfSummary(milexp2[milexp2$T1 == 1,]), method = "render")


## T1 = 0
# Experiment 1
print(dfSummary(milexp1[milexp1$T1 == 0,]), method = "render")

# Experiment 2
print(dfSummary(milexp2[milexp2$T1 == 0,]), method = "render")



### Balance Test (with regard to the main treatment, T1; global model)
milexp_b <- na.omit(milexp)
covs <- subset(milexp_b, select = c(age, female, rank, edu, hometown,
                                  sibling, anger, risk, moral, dove,
                                  partyIDstrength))
bal0 <- bal.tab(T1 ~ covs, data = milexp_b)
love.plot(bal0, binary = "std", thresholds = c(m = .1))
```

(5) Lastly, we present the main analysis results (H1, H2) in a table with stargazer.

```{r, error=FALSE, warning=FALSE, message=FALSE}
H1_1 = lm(Value ~ T1 + Variable + study, data = datstack)
H2_1 = lm(Value ~ T1*T2 + Variable + study, data = datstack)
stargazer(H1_1, H2_1, se = starprep(H1_1, H2_1), type = "latex")
```